Truncated Pareto Distribution (truncpareto) — Bounded Power Laws#
The truncated Pareto is a power-law distribution with an upper cutoff. It’s useful when values are heavy-tailed (many small, a few very large) but a true “infinite tail” is impossible due to physics, policy, or measurement limits.
What you’ll learn#
classification, support, and parameter constraints
the PDF/CDF/PPF in closed form (and why truncation matters)
moments (mean/variance/skewness/kurtosis), characteristic function, and entropy
parameter intuition: how
b(tail exponent) andc(upper cutoff) change the shapea NumPy-only sampler via inverse CDF + Monte Carlo validation
practical usage via
scipy.stats.truncpareto(pdf,cdf,rvs,fit)
import numpy as np
import plotly.express as px
import plotly.graph_objects as go
import os
import plotly.io as pio
from scipy import optimize, stats
# Plotly rendering (CKC convention)
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")
# Reproducibility
rng = np.random.default_rng(7)
np.set_printoptions(precision=4, suppress=True)
1) Title & Classification#
Name:
truncpareto(upper truncated Pareto distribution)Type: Continuous
Support (standardized): \(x \in [1, c]\)
Parameter space (standardized): \(b > 0\), \(c > 1\)
SciPy also supports a shift/scale: if \(Y = \mathrm{loc} + \mathrm{scale}\,X\), then
We write:
2) Intuition & Motivation#
2.1 What it models#
A Pareto/power-law model says “small values are common, large values are rare, and there’s no typical scale.” Real systems often violate the infinite tail assumption: there is usually a maximum feasible value.
truncpareto captures this by taking a Pareto-like density on \([1, \infty)\) and conditioning on \(X \le c\):
2.2 Typical real-world use cases#
Earthquake energy / event sizes: power-law behavior with a physical upper cutoff.
Wildfire sizes: heavy-tailed, but bounded by geography and time.
File sizes / network flows: heavy tails with protocol or system caps.
Wealth / claims with policy caps: heavy-tailed outcomes with maximum payout.
2.3 Relations to other distributions#
Pareto: as \(c \to \infty\), the truncation vanishes and
truncparetoapproaches a standard Pareto.Log-uniform: as \(b \to 0^+\), the density approaches \(f(x) \propto 1/x\) on \([1,c]\) (uniform in \(\log x\)).
Alternatives: lognormal, Weibull, and exponential families can mimic heavy tails over finite ranges; model comparison is often empirical.
3) Formal Definition#
Let \(b>0\) and \(c>1\). Define the normalization constant
3.1 PDF#
3.2 CDF#
3.3 Quantile function (PPF)#
Solving \(F(x)=u\) gives
def _check_truncpareto_params(b: float, c: float) -> tuple[float, float]:
b = float(b)
c = float(c)
if not (b > 0):
raise ValueError("b must be > 0")
if not (c > 1):
raise ValueError("c must be > 1")
return b, c
def truncpareto_logpdf(x: np.ndarray, b: float, c: float) -> np.ndarray:
b, c = _check_truncpareto_params(b, c)
x = np.asarray(x, dtype=float)
logc = np.log(c)
den = -np.expm1(-b * logc) # 1 - c^{-b}
log_k = np.log(b) - np.log(den)
out = np.full_like(x, -np.inf, dtype=float)
mask = (x >= 1) & (x <= c)
out[mask] = log_k - (b + 1) * np.log(x[mask])
return out
def truncpareto_pdf(x: np.ndarray, b: float, c: float) -> np.ndarray:
return np.exp(truncpareto_logpdf(x, b, c))
def truncpareto_cdf(x: np.ndarray, b: float, c: float) -> np.ndarray:
b, c = _check_truncpareto_params(b, c)
x = np.asarray(x, dtype=float)
logc = np.log(c)
den = -np.expm1(-b * logc) # 1 - c^{-b}
out = np.zeros_like(x, dtype=float)
out[x >= c] = 1.0
mask = (x >= 1) & (x < c)
if np.any(mask):
num = -np.expm1(-b * np.log(x[mask])) # 1 - x^{-b}
out[mask] = num / den
return out
def truncpareto_ppf(u: np.ndarray, b: float, c: float) -> np.ndarray:
b, c = _check_truncpareto_params(b, c)
u = np.asarray(u, dtype=float)
if np.any((u < 0) | (u > 1)):
raise ValueError("u must be in [0, 1]")
c_neg_b = np.exp(-b * np.log(c)) # c^{-b}
x_neg_b = 1.0 - u * (1.0 - c_neg_b)
x = x_neg_b ** (-1.0 / b)
return np.clip(x, 1.0, c)
4) Moments & Properties#
Because the support is bounded (\(1\le X\le c\)), all moments exist.
4.1 Raw moments#
For any real \(r\):
This gives the closed form
In particular, the mean is \(\mathbb{E}[X]=\mathbb{E}[X^1]\) and the second raw moment is \(\mathbb{E}[X^2]\). The variance is \(\mathrm{Var}(X)=\mathbb{E}[X^2]-\mathbb{E}[X]^2\).
4.2 Skewness and kurtosis#
Using raw moments \(m_r=\mathbb{E}[X^r]\):
Then
4.3 MGF / characteristic function#
Since \(X\) is bounded, the MGF exists for all real \(t\):
It can be written in terms of the upper incomplete gamma function (for \(t\ne 0\)):
The characteristic function is \(\varphi_X(\omega)=M_X(i\omega)\).
4.4 Entropy#
The differential entropy has a closed form:
(A useful numerically stable identity is \(\dfrac{c^{-b}}{1-c^{-b}} = \dfrac{1}{c^b-1}\).)
def truncpareto_raw_moment(order: float, b: float, c: float) -> float:
b, c = _check_truncpareto_params(b, c)
r = float(order)
logc = np.log(c)
den = -np.expm1(-b * logc) # 1 - c^{-b}
k = b / den
if np.isclose(r, b):
return k * logc
num = np.expm1((r - b) * logc) # c^{r-b} - 1
return k * num / (r - b)
def truncpareto_moments(b: float, c: float) -> dict:
m1 = truncpareto_raw_moment(1.0, b, c)
m2 = truncpareto_raw_moment(2.0, b, c)
m3 = truncpareto_raw_moment(3.0, b, c)
m4 = truncpareto_raw_moment(4.0, b, c)
var = m2 - m1**2
std = np.sqrt(var)
mu3 = m3 - 3 * m1 * m2 + 2 * m1**3
mu4 = m4 - 4 * m1 * m3 + 6 * (m1**2) * m2 - 3 * m1**4
skew = mu3 / (std**3)
excess_kurt = mu4 / (var**2) - 3
return {
"mean": m1,
"var": var,
"skew": skew,
"excess_kurt": excess_kurt,
"mode": 1.0,
"median": float(truncpareto_ppf(0.5, b, c)),
}
def truncpareto_entropy(b: float, c: float) -> float:
b, c = _check_truncpareto_params(b, c)
logc = np.log(c)
den = -np.expm1(-b * logc) # 1 - c^{-b}
log_k = np.log(b) - np.log(den)
# E[log X] = 1/b - log(c)/(c^b - 1)
elogx = (1.0 / b) - (logc / np.expm1(b * logc))
return -log_k + (b + 1) * elogx
# Quick numerical check against SciPy
b0, c0 = 2.0, 5.0
ours = truncpareto_moments(b0, c0)
mean_s, var_s, skew_s, exkurt_s = stats.truncpareto.stats(b0, c0, moments="mvsk")
{
"mean (ours, scipy)": (ours["mean"], float(mean_s)),
"var (ours, scipy)": (ours["var"], float(var_s)),
"skew (ours, scipy)": (ours["skew"], float(skew_s)),
"excess_kurt (ours, scipy)": (ours["excess_kurt"], float(exkurt_s)),
"entropy (ours, scipy)": (truncpareto_entropy(b0, c0), float(stats.truncpareto.entropy(b0, c0))),
}
{'mean (ours, scipy)': (1.666666666666667, 1.6666666666666667),
'var (ours, scipy)': (0.5752178731265971, 0.5752178731265976),
'skew (ours, scipy)': (1.897052928779272, 1.8970529287792715),
'excess_kurt (ours, scipy)': (3.587240444439618, 3.5872404444395976),
'entropy (ours, scipy)': (0.5648510858655371, 0.564851085865537)}
5) Parameter Interpretation#
truncpareto has two standardized parameters:
\(b\) (shape / tail exponent):
Appears in \(x^{-(b+1)}\).
Larger \(b\) puts more mass near 1 (lighter tail).
Smaller \(b\) produces a heavier tail (more probability near \(c\)).
Limit: \(b \to 0^+\) gives \(f(x) \to \dfrac{1}{x\log c}\) (log-uniform on \([1,c]\)).
\(c\) (upper cutoff):
The maximum possible value (standardized).
Larger \(c\) increases the range and typically increases mean/variance.
Limit: \(c\to\infty\) recovers the classic Pareto on \([1,\infty)\).
Scaling to arbitrary bounds#
If you want support \([x_{\min}, x_{\max}]\) with \(0<x_{\min}<x_{\max}\), set
In SciPy this corresponds to scale=x_min, c=x_max/x_min (and typically loc=0).
# How b changes the shape (fixed c)
c_fixed = 10.0
x = np.linspace(1.0, c_fixed, 800)
b_values = [0.3, 1.0, 2.0, 6.0]
fig = go.Figure()
for b in b_values:
fig.add_trace(
go.Scatter(x=x, y=truncpareto_pdf(x, b, c_fixed), mode="lines", name=f"b={b:g}")
)
fig.update_layout(
title=f"truncpareto PDF: varying b (c={c_fixed:g})",
xaxis_title="x",
yaxis_title="density",
width=900,
height=420,
)
fig.show()
# How c changes the shape (fixed b)
b_fixed = 2.0
c_values = [1.5, 2.0, 5.0, 20.0]
x = np.linspace(1.0, max(c_values), 900)
fig = go.Figure()
for c in c_values:
mask = x <= c
fig.add_trace(
go.Scatter(
x=x[mask],
y=truncpareto_pdf(x[mask], b_fixed, c),
mode="lines",
name=f"c={c:g}",
)
)
fig.update_layout(
title=f"truncpareto PDF: varying c (b={b_fixed:g})",
xaxis_title="x",
yaxis_title="density",
width=900,
height=420,
)
fig.show()
6) Derivations#
6.1 Normalization constant#
Start with \(f(x)=Kx^{-(b+1)}\) on \([1,c]\). Enforce total probability 1:
So
6.2 Expectation (general power moment)#
For \(r\in\mathbb{R}\),
If \(r\ne b\):
If \(r=b\), the integrand becomes \(x^{-1}\) and the integral is \(\log c\).
6.3 Variance#
Using the raw moments \(m_1=\mathbb{E}[X]\) and \(m_2=\mathbb{E}[X^2]\):
6.4 Likelihood (i.i.d. sample)#
For data \(x_1,\dots,x_n\) with \(1\le x_i\le c\):
The log-likelihood is
The constraint \(c \ge \max_i x_i\) is part of the model’s support. Because the only dependence on \(c\) is through \(-n\log(1-c^{-b})\) (which decreases as \(c\) increases), the MLE is
With \(c\) fixed/known, the score equation for \(b\) is
which is typically solved numerically.
def truncpareto_loglikelihood(x: np.ndarray, b: float, c: float) -> float:
x = np.asarray(x, dtype=float)
b, c = _check_truncpareto_params(b, c)
if np.any((x < 1) | (x > c)):
return -np.inf
n = x.size
logc = np.log(c)
den = -np.expm1(-b * logc) # 1 - c^{-b}
log_k = np.log(b) - np.log(den)
return n * log_k - (b + 1) * np.sum(np.log(x))
def truncpareto_mle_c(x: np.ndarray) -> float:
x = np.asarray(x, dtype=float)
if np.any(x < 1):
raise ValueError("For the standardized model, all data must satisfy x >= 1.")
return float(np.max(x))
def _score_b(b: float, x: np.ndarray, c: float) -> float:
n = x.size
sum_logx = np.sum(np.log(x))
return (n / b) - (n * np.log(c) / (c**b - 1.0)) - sum_logx
def truncpareto_mle_b(x: np.ndarray, c: float) -> float:
x = np.asarray(x, dtype=float)
c = float(c)
if not (c > 1):
raise ValueError("c must be > 1")
if np.any((x < 1) | (x > c)):
raise ValueError("All data must lie in [1, c] for the standardized model.")
lo = 1e-8
hi = 10.0
f_lo = _score_b(lo, x, c)
f_hi = _score_b(hi, x, c)
# If the score is already negative at tiny b, the likelihood is decreasing from the boundary.
if f_lo <= 0:
return lo
# Expand the upper bracket until the score becomes negative.
while (f_hi > 0) and (hi < 1e6):
hi *= 2
f_hi = _score_b(hi, x, c)
if f_hi > 0:
# Extremely flat / boundary-ish case
return hi
return float(optimize.brentq(lambda bb: _score_b(bb, x, c), lo, hi))
# Demo: MLE on synthetic data
b_true, c_true = 2.0, 5.0
x = stats.truncpareto.rvs(b_true, c_true, size=4000, random_state=rng)
c_hat = truncpareto_mle_c(x) # boundary MLE
b_hat_known_c = truncpareto_mle_b(x, c_true)
{
"true (b, c)": (b_true, c_true),
"MLE c_hat (=max x)": c_hat,
"MLE b_hat (assuming c known)": b_hat_known_c,
}
{'true (b, c)': (2.0, 5.0),
'MLE c_hat (=max x)': 4.966201346197002,
'MLE b_hat (assuming c known)': 1e-08}
7) Sampling & Simulation#
The CDF has a closed-form inverse, so inverse transform sampling is straightforward.
Starting from
set \(U\sim\mathrm{Uniform}(0,1)\) and solve \(F(X)=U\):
This is efficient, stable, and uses only uniform random numbers.
def truncpareto_rvs_numpy(b: float, c: float, size: int, rng: np.random.Generator) -> np.ndarray:
b, c = _check_truncpareto_params(b, c)
u = rng.random(size)
return truncpareto_ppf(u, b, c)
# Compare NumPy-only sampler to SciPy sampler (quick two-sample KS test)
b0, c0 = 2.0, 5.0
n = 50_000
samples_numpy = truncpareto_rvs_numpy(b0, c0, size=n, rng=rng)
samples_scipy = stats.truncpareto.rvs(b0, c0, size=n, random_state=rng)
stats.ks_2samp(samples_numpy, samples_scipy)
KstestResult(statistic=0.006680000000000019, pvalue=0.21359969753360408, statistic_location=1.5599207153635837, statistic_sign=-1)
8) Visualization#
We’ll visualize:
the PDF and CDF for chosen parameters
Monte Carlo samples (histogram + theoretical PDF overlay)
the empirical vs theoretical CDF
b0, c0 = 1.7, 8.0
x = np.linspace(1.0, c0, 900)
# PDF and CDF
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=truncpareto_pdf(x, b0, c0), mode="lines", name="PDF"))
fig.update_layout(
title=f"truncpareto PDF (b={b0:g}, c={c0:g})",
xaxis_title="x",
yaxis_title="density",
width=900,
height=380,
)
fig.show()
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=truncpareto_cdf(x, b0, c0), mode="lines", name="CDF"))
fig.update_layout(
title=f"truncpareto CDF (b={b0:g}, c={c0:g})",
xaxis_title="x",
yaxis_title="F(x)",
width=900,
height=380,
)
fig.show()
# Monte Carlo sample vs PDF
n = 40_000
samples = truncpareto_rvs_numpy(b0, c0, size=n, rng=rng)
fig = go.Figure()
fig.add_trace(
go.Histogram(
x=samples,
histnorm="probability density",
nbinsx=60,
name="samples",
opacity=0.35,
)
)
fig.add_trace(go.Scatter(x=x, y=truncpareto_pdf(x, b0, c0), mode="lines", name="theoretical PDF"))
fig.update_layout(
title=f"Monte Carlo samples vs PDF (n={n:,}, b={b0:g}, c={c0:g})",
xaxis_title="x",
yaxis_title="density",
width=900,
height=420,
bargap=0.02,
)
fig.show()
# Empirical CDF vs theoretical CDF
xs = np.sort(samples)
ecdf = np.arange(1, xs.size + 1) / xs.size
fig = go.Figure()
fig.add_trace(go.Scatter(x=xs, y=ecdf, mode="lines", name="empirical CDF"))
fig.add_trace(go.Scatter(x=x, y=truncpareto_cdf(x, b0, c0), mode="lines", name="theoretical CDF"))
fig.update_layout(
title="Empirical vs theoretical CDF",
xaxis_title="x",
yaxis_title="F(x)",
width=900,
height=420,
)
fig.show()
9) SciPy Integration (scipy.stats.truncpareto)#
SciPy implements the standardized distribution with shape parameters b and c.
Key methods:
pdf,logpdf,cdf,ppfrvs(sampling)stats(..., moments='mvsk')andentropyfit(MLE)
Remember: with loc and scale, SciPy transforms \(X\) via \(Y=\mathrm{loc}+\mathrm{scale}\,X\).
b0, c0 = 1.7, 8.0
rv = stats.truncpareto(b0, c0)
# Evaluate core functions
x = np.linspace(1.0, c0, 6)
{
"x": x,
"pdf": rv.pdf(x),
"cdf": rv.cdf(x),
"ppf(0.1,0.5,0.9)": rv.ppf([0.1, 0.5, 0.9]),
}
# Sampling
data = rv.rvs(size=5000, random_state=rng)
# Fitting (fix loc=0, scale=1 for the standardized model)
b_hat, c_hat, loc_hat, scale_hat = stats.truncpareto.fit(data, floc=0, fscale=1)
(b_hat, c_hat, loc_hat, scale_hat)
(1.6649106936679934, 7.94363698958849, 0, 1)
10) Statistical Use Cases#
10.1 Hypothesis testing / goodness-of-fit#
A common workflow is:
fit \((b,c)\) (often fixing
loc/scale)assess fit using diagnostic plots (PDF/CDF/QQ) and tests (e.g. KS)
Important: when parameters are estimated from the same data, classical p-values from KS tests are no longer exact. A more principled approach is a parametric bootstrap under the fitted model.
10.2 Bayesian modeling#
Truncated Pareto priors are natural when you want a heavy-tailed prior but also know a hard upper bound. A classic example is an unknown “maximum” parameter \(\theta\) (e.g., the endpoint of a Uniform distribution).
10.3 Generative modeling#
If you need realistic synthetic data with a capped power-law tail (event sizes, file sizes, claim sizes), truncpareto is a simple, interpretable generator.
# 10.1 Example: (diagnostic) KS test against a fitted truncpareto
b_true, c_true = 1.8, 10.0
x = stats.truncpareto.rvs(b_true, c_true, size=2500, random_state=rng)
b_hat, c_hat, _, _ = stats.truncpareto.fit(x, floc=0, fscale=1)
ks = stats.kstest(x, "truncpareto", args=(b_hat, c_hat))
{
"true (b, c)": (b_true, c_true),
"fitted (b, c)": (float(b_hat), float(c_hat)),
"KS statistic": float(ks.statistic),
"KS p-value (not exact after fit)": float(ks.pvalue),
}
{'true (b, c)': (1.8, 10.0),
'fitted (b, c)': (1.829190701746747, 9.973639795656418),
'KS statistic': 0.012295004694353129,
'KS p-value (not exact after fit)': 0.8395284823777183}
# 10.2 Example: truncated Pareto prior for an unknown Uniform endpoint
# Data model: Y_i | theta ~ Uniform(0, theta)
# Prior: theta ~ TruncPareto(b0, c=U/L) scaled to [L, U]
theta_true = 8.0
n = 40
y = rng.uniform(0.0, theta_true, size=n)
y_max = float(np.max(y))
L, U = 1.0, 20.0
b_prior = 2.0
prior = stats.truncpareto(b_prior, U / L, loc=0.0, scale=L)
# Posterior under the Uniform likelihood:
# p(theta | y) ∝ theta^{-n} * theta^{-(b_prior+1)} = theta^{-(b_prior+n+1)} on [max(L, y_max), U]
L_post = max(L, y_max)
b_post = b_prior + n
post = stats.truncpareto(b_post, U / L_post, loc=0.0, scale=L_post)
grid = np.linspace(L, U, 900)
fig = go.Figure()
fig.add_trace(go.Scatter(x=grid, y=prior.pdf(grid), mode="lines", name=f"prior (b={b_prior:g})"))
fig.add_trace(go.Scatter(x=grid, y=post.pdf(grid), mode="lines", name=f"posterior (b={b_post:g})", line=dict(width=3)))
fig.add_vline(x=y_max, line_dash="dash", line_color="gray", annotation_text="max(y)")
fig.update_layout(
title=f"Uniform endpoint inference (n={n}, theta_true={theta_true:g})",
xaxis_title="theta",
yaxis_title="density",
width=900,
height=420,
)
fig.show()
{
"y_max": y_max,
"posterior mean": float(post.mean()),
"95% credible interval": tuple(map(float, post.ppf([0.025, 0.975]))),
}
{'y_max': 7.96832887580269,
'posterior mean': 8.162678360578365,
'95% credible interval': (7.9731336719402, 8.699845415743587)}
# 10.3 Example: generative modeling + CCDF (survival) on log-log axes
b0, c0 = 1.3, 50.0
n = 80_000
s = truncpareto_rvs_numpy(b0, c0, size=n, rng=rng)
s_sorted = np.sort(s)
# Empirical survival S(x) = P(X > x)
surv_emp = 1.0 - (np.arange(1, n + 1) / n)
x_grid = np.logspace(0, np.log10(c0), 500)
surv_theory = 1.0 - truncpareto_cdf(x_grid, b0, c0)
fig = go.Figure()
fig.add_trace(go.Scatter(x=s_sorted, y=surv_emp, mode="lines", name="empirical survival"))
fig.add_trace(go.Scatter(x=x_grid, y=surv_theory, mode="lines", name="theoretical survival"))
fig.update_layout(
title=f"Survival function on log-log scale (b={b0:g}, c={c0:g})",
xaxis_title="x",
yaxis_title="P(X > x)",
xaxis_type="log",
yaxis_type="log",
width=900,
height=420,
)
fig.show()
11) Pitfalls#
Invalid parameters: require \(b>0\) and \(c>1\).
Support mismatches: the standardized model assumes \(x\in[1,c]\). If your natural lower bound is \(x_{\min}\ne 1\), use a scale transform.
Boundary behavior when fitting:
With \(c\) unknown, the MLE is the sample maximum (\(\hat c=\max x_i\)), which is a boundary estimate.
If you know the physical cutoff, treat \(c\) as fixed rather than estimated.
Numerical issues when \(c\approx 1\): \(1-c^{-b}\) can be tiny; prefer
logpdfandexpm1/log1p-style computations.Interpreting power laws: many distributions can look linear on log-log plots over a short range; validate with diagnostics and model comparison.
12) Summary#
truncparetois a continuous bounded power-law on \([1,c]\) with parameters \(b>0\) and \(c>1\).PDF: \(f(x)=\dfrac{b}{1-c^{-b}}x^{-(b+1)}\) on \([1,c]\); CDF: \(F(x)=\dfrac{1-x^{-b}}{1-c^{-b}}\).
All moments exist (bounded support); raw moments have a simple closed form.
Sampling is easy via the inverse CDF: \(X=(1-U(1-c^{-b}))^{-1/b}\).
SciPy’s
scipy.stats.truncparetoprovidespdf/cdf/rvs/fit, plus moments and entropy.